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Abstract 


1.0 Introduction 


An efficient incremental -iterative approach for dif- 
ferentiating advanced flow codes is successfully demon- 
strated on a 2D inviscid model problem. The method 
employs the reverse-mode capability of the automatic- 
differentiation software tool ADI FOR 3.0 , and is 
proven to yield accurate first-order aerodynamic sensi- 
tivity derivatives. A substantial reduction in CPU time 
and computer memory is demonstrated in comparison 
with results from a straight-forward, black-box reverse - 
mode application of ADI FOR 3.0 to the same flow 
code. An ADIFOR-assisted procedure for accurate sec- 
ond-order aerodynamic sensitivity derivatives is suc- 
cessfully verified on an inviscid transonic lifting airfoil 
example problem. The method requires that first-order 
derivatives are calculated first using both the forward 
(direct) and reverse (adjoint) procedures; then, a very 
efficient non-iterative calculation of all second-order 
derivatives can be accomplished . Accurate second de- 
rivatives ( i. e. , the complete Hessian matrices) of lift , 
wave-drag , and pitching-moment coefficients are calcu- 
lated with respect to geometric-shape , angle-of-attack , 
and freestream Mach number 


Computing sensitivity derivatives (SDs) from high- 
fidelity, nonlinear CFD codes is an enabling technology 
for design of advanced concept vehicles. In recent years 
significant progress has been achieved in the efficient 
calculation of accurate SDs from these CFD codes 1 . 
The automatic differentiation (AD) software tool 
ADIFOR (Automatic Differentiation of FORTRAN) 
has been proven an effective tool for extracting aerody- 
namic SDs from these modern CFD codes 2 ' 6 . The pre- 
sent study will essentially build on earlier studies 3 ' 6 in 
an effort to exploit the full potential of the latest version 
of ADIFOR 3.0 7 for obtaining SDs from CFD codes. 

In Ref. 2, a strategy was first proposed and later 
successfully demonstrated in Ref. 3, whereby AD was 
applied to a CFD code in incremental-iterative (I-I) 
form. This hybrid scheme (known as the ADII method) 
was designed to achieve in part the computational effi- 
ciency of a hand-differentiated (HD) approach, and at 
the same time capture identical accuracy while main- 
taining (at least in part) the ease-of-implementation of a 
straightforward black-box (BB) application of AD. The 
2D effort of Ref. 3 was later extended to the 3D code 
CFL3D, including “in-parallel” computation of the de- 
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rivatives 8,9 . Appropriate references to the version of 
CFL3D used can be found in Ref. 6. 

The success reported in these previous works 3,8,9 
could be considered limited, however, because all 
ADIFOR implementations reported therein were “for- 
ward-mode” (direct) differentiations. It is very difficult 
to make any forward-mode implementation of deriva- 
tive calculations computationally competitive with a 
“reverse-mode” (adjoint) implementation whenever the 
number of design variables (NDV) of interest is consid- 
erably larger than the number of output functions 
(NOF) of interest; and NDV much greater than NOF is 
more typical with aerodynamic design problems. In 
recent studies, the new reverse-mode capability of 
ADIFOR 3.0 (not available for the earlier referenced 
studies) has been successfully verified by application 0 
to an “in-parallel” version of CFL3D and application 10 
to a sequential linear aerodynamics code, resulting in 
accurate design SDs as well as stability and control 
derivatives, respectively. The application reported in 
Ref. 6 involved BB AD of the entire CFD code, but 
iterative execution of the reverse-mode was over only 
the last function iteration tree. 

In the present study, it is proposed and demon- 
strated that the reverse-mode capability of ADIFOR 3.0 
can also be applied to CFD codes in I-I form, resulting 
in a hybrid adjoint-variable (AV) scheme (known 
herein as the ADII-AV method) that is analogous to the 
forward-mode ADII scheme of Ref. 3 and elsewhere. 
The motivation of this new reverse-mode ADII-AV 
scheme is identical to that of the earlier forward-mode 
ADII method: greater computational efficiency is 
sought over a BB implementation of AD, without any 
loss of accuracy in the calculated SDs and without un- 
manageable complications upon implementation. 

Following development of the proposed new ADII- 
AV scheme, the second focus of the present study is 
that of calculating second-order (SO) aerodynamic SDs 
from CFD codes. This second part of the present study 
is another extension of Ref. 3, wherein the computa- 
tional issues associated with calculating these higher- 
order derivatives were addressed, and sample calcula- 
tions of SO derivatives using AD were reported from a 
2D CFD code. In Ref. 3, four procedures for calculat- 
ing SO CFD SDs were proposed, but only one of the 
less efficient methods was actually tested; it should also 
be noted here that ADIFOR 3.0 currently provides three 
forward-mode variations for the calculation of SO SDs 
by similarly inefficient methods. The most efficient (for 
large NDV) SO SD scheme was not tested in the earlier 
study 3 , but has been successfully implemented in the 
present study. Reverse-mode adjoint-based differentia- 
tion is required within this efficient SO SD scheme; 
therefore, with the availability of ADIFOR 3.0 and the 


new ADII-AV scheme, the door has been opened for 
implementation and testing of the most efficient SO SD 
scheme for CFD codes. The results of this effort to date 
are reported herein. In a companion study 11 these effi- 
ciently computed SO SDs have been used to demon- 
strate an approach for CFD input uncertainty propaga- 
tion and robust design optimization for a quasi- ID flow 
application. 

2.0 Basic Equations and 
Theoretical Development 

The equations summarized subsequently are dis- 
cussed in greater detail in the references, in particular, 
Ref. 3. These concepts are known in the mathematical 
optimization community (Ref. 12) but the details devel- 
oped here do not appear to be generally known 
throughout the CFD community. The aerodynamic 
output functions of interest, F, and the discretized con- 
servation laws of steady compressible fluid flow, R , 
including boundary conditions can be represented sym- 
bolically as 

F =F{Q{b),X{b\b) (1) 

(aerodynamic output functions) 


R = R{Q{b),X{b),b) = 0 (2) 

(nonlinear stale equations) 

where Q is the vector of state (Field) variables, X is the 
vector of computational grid coordinates, and b is the 
vector of input (design) variables. 


2.1 First-Order Sensitivity Derivatives 


By the following preliminary definitions, index 
(summation) notation is now' introduced. (This notation 
will be necessary to avoid subsequent ambiguity when 
the SO SD methods are presented.) 


y db . 


f dF^ 

db 


R'„ = 


dR , ( dR 


J'i 


db . 




()’ _ d Q,n 

^ db. 


r dQ) 

db 


, dX„ 

y f p 

M ^ 


db 


v 


Jpj 


The forward-mode (direct) approach for calculating 
first-order (FO) SDs is developed by differentiation of 
Eq. (1) and (2) with respect to the design variables; the 
result is 


’ db, M K 




A 

dX. 


df 


PJ+ db, 


(3) 
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, _ dR, _ dR, 


aa. 


+ M 


dR, , dR, 

x « + ^ =0 « (i) 


In the above, i, 7 , and / are “free” indices, and repeated 
indices m and p are (by convention) “summation” indi- 
ces. The reverse-mode (adjoint) approach for the FO 
SDs is developed in a non-conventional manner starting 
with application of the chain rule 

dF, dR, dF t 

1 ~ = (5) 

dR, dQ m dQ m 

With Eq. (5), it then follows from eq. (4) that 


where the superscript N is the iteration (pseudo-time 
step) index, and the operator 


pN _ _ 

r ml ~~ - ~ ” 


f 21 n 'V \ 


a/?; v 


dR f 

da, 


(id 


represents the solution algorithm of the particular CFD 
code of choice. The tilde (~) in Eq. (II) serves to indi- 
cate that P^ ( can be viewed as any computationally 


efficient approximation (often a very crude approxima- 
tion) of the exact operator associated with true Newton- 
Raphson iteration. Thus the CFD solution algorithm is 
simply quasi-Newton iteration. 


dpi ( y _ dp dR , , 

dQ„ mi dR,dQ m ^ 


dR, 


f 


dR, 


dR, 


Y' + 

dx„ Xpi db 


J 


( 6 ) 


and Eq. (3) becomes 




+ 


( d R i Y ' + ^ 

dT PJ db. 


_dp 

dR, 

dp_ 

dx Xpi db, 




(7) 


Comparison of Eq. (5) and (7) with a more conven- 
tional development of the reverse-mode method yields 
the identity 


.^ = 4 , 

dR, " 


( 8 ) 


where A if is more commonly called the “adjoint vari- 
able.” Using Eq. ( 8 ), the more conventional presenta- 
tion of Eq. (5) is 


The I-I method for solving the forward-mode, FO 
SD Eq. (4) is 


G 'M +1 xVM — p 

mj r ml n lj 


02 ) 


where superscript M is the FO SD iteration index, and 


K 


M 


= ^Plo' m +^-X' i dfr 
dQ„ mj dx p pj db, 


(13) 


With the I-I methodology, the CFD flow solution op- 
erator P ml is also used to solve the SD equations; this 

operator in Eq. (12) is evaluated and fixed using the 
steady-state solution for the nonlinear flow. The requi- 
site terms of Eq. (13) are constructed either by hand 
differentiation (i.e., the HDTI method, which is very 
tedious and time consuming to complete with accuracy 
for advanced CFD codes) or by AD, which is the for- 
ward-mode ADII method of previous studies. 

In contrast with the ADII method, a straightforward 
BB application of AD to the CFD code, which is the 
ADBB method, is represented symbolically as 

mj =Qmj - P ml R ,j ~ P n tl , R l ( 1 4 ) 


G , m = K 


dR, dF, 
dQ„, dQ m 


(9) 


One objective of this particular development of the AV 
method for aerodynamic SDs is to ensure that the rela- 
tionship given by Eq. ( 8 ) is clearly understood. 

The I-I strategies for solving the preceding equa- 
tions are reviewed here; additional detail is found in 
Ref 3 and elsewhere. The I-I method for solving the 
nonlinear flow of Eq. (2) is 


Q n+] = Q' 

m x-' n 


n N n N 

Fnl R l 


( 10 ) 


Clearly ADII (Eq. (12) and (13)) and ADBB (Eq. (14)) 
yield the same result at steady-state convergence of 
each (recall Eq. (2)); however, ADII is potentially more 
efficient than ADBB due to user intervention in the 
application of AD. With ADII 

1. The operator P^ f can be evaluated only once using 

the steady-state field variables Q and then reused 
for all M iterations and for all j = NDV design 

variables in obtaining the O' . 

2 . All derivatives except Q' nj can be computed once 

outside the iteration loop and frozen for reuse in- 
side the loop. 
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3. Evaluation of the terms in Eq. (!4) can he 

avoided completely, for all iterations and all design 
variables. 

The I-I method for solving the reverse-mode, AV, 
FO SD Eq. (9) is 


required for efficient evaluation of the term 


iW 


dR, 

d<2„, 


inEq. (16). 


2.2 Second-Order Sensitivity Derivatives 


where 


+1 l\f P 

^il ~ A'H ~ 


=^l— l - + ^L 

30 . 30 . 


(15) 


(16) 


The requisite terms of Eq. (16) are constructed ei- 
ther by hand (i.e., the HDII-AV method, having the 
same drawbacks as the forward- mode HDIT method), or 
by AD, which is the proposed new ADII-AV scheme. 
The BB AD in reverse-mode (the ADBB-AV method) 
has been verified in Ref. 6. The objective of the pro- 
posed ADII-AV scheme is improved computational 
efficiency over the ADBB-AV approach without result- 
ing loss of accuracy or significant loss in the ease-of- 
implementation. The mechanisms from which improved 
computational efficiency can be expected are analogous 
to that explained previously when the forward-mode 
ADII and ADBB methods were contrasted. Further- 
more, the ADII-AV scheme should lend itself to more 
permanent generalized coding implementations than the 
ADBB-AV approach. This is because with the 
ADII-AV method, the manner in which AD is applied 
is independent of, yet valid for, all the particular aero- 
dynamic inputs and outputs of interest. 

The forward -mode application of ADIFOR has al- 
ways been very effective in constructing FORTRAN 
source code for the computationally efficient, repeated 
calculation of the vector (or matrix) product that results 
from the post- multiplication of a large Jacobian matrix 
by a known input vector (or matrix). This attribute of 
forward-mode AD is exactly what was required to con- 
struct the ADII method; specifically, the terms 

— Q M and X ni of Eq. (13) are of this type. 

30. ' asr, " 

In contrast, however, the forward-mode application of 
ADIFOR constructs source code which is prohibitively 
inefficient for calculating the pre-multiplication of a 
large Jacobian matrix by a known input vector (or ma- 
trix). This weakness of the forward-mode application of 
ADIFOR is exactly the strength of the reverse-mode 
option now available in ADIFOR 3.0. Thus, the pro- 
posed new efficient ADII-AV scheme has become pos- 
sible with this reverse-mode capability. That is, through 
reverse-mode application of ADIFOR 3.0, it is now 
feasible to construct (automatically) the source code 


The SO SD methods are presented subsequently us- 
ing the index notation and beginning with the following 
preliminary definitions: 



db k db j 


( d 2 F 


\ 



Jijk 




d% 

db.db, 


db k dbj 


(d 2 R ^ 

V 


db 2 


jtjk 


(djOj 

db 2 


Jmjk 


y ~ 

A PJk ” 


£L_ 

db k dbj 


fd 2 X^ 


db 2 


Jpjk 


v _ d\ f _ 

4/* — „ ~ 

d h 


f— 1 

ydb) !k 


in 

db k dRf 


GL = 


dG „ 
db t 


f dG_ S 
db 


Jimk 


The following differential operator is also introduced 
for subsequent notational compactness: 


Db k dQ, Qnk dx q * db k 


(17) 


where repeated indices n and q are summation indices. 

Differentiation of the FO forward-mode Eq. (3) 
and (4) with respect to the design variables yields SO 
Method 1 


f*=jLL. 

,}k db.db 


= M-Q’ 

• .V + Bb 


( 18 ) 


sx. 


k pjk 


Db, 
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D” _ d R, _ c)R t 

lX ljk - „ „ ~ V£| 


dh t dbj 30 , 

DR' 


B# + ax“ 




/y* 


(19) 


+ ■ 


Db, 


■ = 0 


m 


The terms of DF' / Db k and DR'. / Db k are many 

and very complicated; detailed expansion of these terms 
is provided in the appendix. Using symmetry of the 

Hessian (2*,^ = Qmkj % SO Method 1 requires 

(NDV NDV)/ 2 solutions of the large linear 

systems of Eq. (19) for Q*,j k ; in addition, the method 

requires NDV solutions of Eq. (4) for the FO SDs Q ' }/ . 

SO Method 1 was verified for a 2D CFD code in Ref. 3 
by ADBB differentiation of the code's existing HDII 
scheme (Eq. ( 1 2) and ( 1 3)) for the FO SDs. 

Alternatively, differentiation of the FO reverse- 
mode, Eq. (7), (8), and (9), with respect to the design 
variables yields SO Method 2 




y* 


db k dbj 

f 


‘ilk 


a/?, , a/?, 

x P j + 


v 


dF i dR, 

a3T "aF 


V*' 
A 


db jj 

DF'; 


J 




o = 


dOj. _ ,, dR, . DO,, 

db « 3 Q. Db k 


= 0 


imk 


( 20 ) 


( 21 ) 


SO Method 2 requires NDV X NOF solutions of the 
large linear systems of Eq. (21) for A' lk ; in addition, 
the method requires NDV solutions of the FO Eq. (4) 
for Q' m j plus NOF solutions of the FO Eq. (9) for X i} . 

This SO Method 2 is eliminated from further considera- 
tion because it is unconditionally less computationally 
efficient than the remaining two SO SD methods. 

Introduction of the AV approach within SO Method 
1 to eliminate Q* ljk yields SO Method 3 


SO Method 4 is similar and computationally equivalent 
to SO Method 3, and is developed by introduction of 
the AV approach within SO Method 2 to eliminate 

A' (k ; the result is the identity 


, D K _ r , DC.. 


Db k 


-mi 


Db K 


(23) 


where SO Method 4 uses Eq. (23) to replace equivalent 
terms within SO Method 3. It is important to note that 
the equivalent SO SD Methods 3 and 4 do not require 
solution of large systems of linear equations for higher- 

order derivatives such as Q*, /k or X ' (k . These two SO 
SD schemes do however require solution of both for- 
ward-mode and reverse-mode Eq. (4) and (9) for Q' m 


"V 


and A i( , respectively. This is a total of only NDV + 
NOF solutions of large systems of linear equations. 


One significant conclusion of the preceding analysis 
is that SO Method 3 or 4 should be computationally 

more efficient whenever NDV " + NDV is greater 
than 2xNOF . With typical design problems in aero- 
dynamics, NDV is often very much larger than NOF, 
and the advantage in favor of Method 3 or 4 for SO SDs 
is then all the greater. Once both the forward-mode and 
reverse-mode schemes are in place for calculating the 
FO SDs, then complete SO SD information is available 
almost “for free”; i.e., the SO SD are obtained through 
an explicit, non-iterative calculation. The source code 
for implementation of Methods 3 or 4 is constructed 
“automatically” via BB application of the forward- 
mode capability of ADIFOR to appropriate pieces of 
the existing source code from which the FO SDs are 
obtained. For example, the extremely complex terms 

DF' / Db . , DR' / Db. , and/or DG. m / Db. (see 

ij k If A. itn k v 

appendix) of Methods 3 and/or 4 are easily constructed 
with a forward-mode application of AD. 

In Ref. 3, SO Methods 3 and 4 were proposed but 
not actually tested. Consequently, one primary goal of 
the present study is successful implementation and veri- 
fication of the highly efficient SO Method 3 (or equiva- 
lently, Method 4); Method 3 is actually chosen in this 
study. 


3,0 Results 


F* = 


^3 F, 3 R, A 

W +l "dx 


P 

DF' 




A,. 


n 

DR' 


Db k Db k 


3.1 First-Order Sensitivity Derivatives, 
ADII-AV Method. Model Problem 

The proposed ADII-AV method (Eq. (15) and (16)) 
has been successfully implemented in a CFD code and 
verified for accuracy on a 2D inviscid internal flow 
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problem. This CFD code solves the 2D Euler equations 
by a conventional upwind finite-volume approach on a 
very coarse grid; but, one that is sufficient for verifying 
SDs. As expected, when FO SDs computed by the new 
ADII-AV scheme are compared with SDs computed by 
a hand-differentiated implementation of Eq. (15) and 
(16) (i.e., the HDII-AV approach), the results are the 
same, iteration-by-iteration. In addition, the accuracy of 
the computed SDs has been successfully verified by a 
finite-difference method. 

Preliminary timings were conducted on a Sun work- 
station to evaluate the potential for improved computa- 
tional efficiency of the new ADII-AV scheme with re- 
spect to the ADBB-AV approach of Ref. 6. Computa- 
tional timing comparisons are given in Table 1 . Rela- 
tive timings are given as CPU time per iteration per 
grid-point per differentiated-aerodynamic-output- 
function. Furthermore, each timing result has been 
scaled by the comparable timing result obtained from 
the very efficient hand-differentiated reverse-mode 
scheme (i.e., the HDII-AV method). 


Table 1. Relative CPU Timing Comparison 


Reverse-Mode Method 
Tested 

Relative Timing* 

ADII-AV / HDII-AV 

4.7 

ADBB-AV /HDII-AV 

7.9 

*CPU Tirne/iteration/grid-poinl/output-function | 


Table 1 illustrates that, although the new ADII-AV 
scheme is almost five times slower than the efficient 
HDII-AV scheme, it represents a substantial improve- 
ment over results obtained from the straight-forward 
black-box procedure (lev; ADBB-AV is about eight 
times slower than HDII-AV). 

Another important computational concern mitigated 
by the new ADII-AV method is computer memory - 
particularly the issue of large disk files created during 
execution of reverse-mode derivative code created by 
ADIFOR 3.0. With the black-box (ADBB-AV) ap- 
proach, these large ADIFOR “log” files (which are cre- 
ated on a forward-pass execution and are read during 
the reverse pass) will accumulate and become larger 
with every iteration of the ADIFOR-enhanced flow 
code. This file growth can rapidly deplete the available 
disk space, even on the largest computers. In Ref. 6 
this difficulty was addressed by development of the 
“iterated reverse-mode” scheme, where only the log 
files for the final forward-pass iteration are stored and 
used during the subsequent iterative solution for the 


derivatives. With the ADII-AV approach, however, the 
required disk space is not as restrictive an issue because 
it remains fixed and does not accumulate during the 
iterative solution process. In the present example, the 
total storage requirement for log files with the ADII- 
AV method is only 64 percent of that required for a 
single iteration of the ADBB-AV method. 

3.2 Second-Order Sensitivity Derivatives, 
SO Method 3, Airfoil Example 

Results are presented subsequently from the suc- 
cessful verification of the proposed efficient non- 
iterative SO Method 3 (Eq. (22)) for computing SO 
SDs. The example problem is steady transonic inviscid 
flow over a NACA 0012 airfoil with freestream Mach 
number (Mj 0.80 and angle-of-attack (a) 1.0 degree. 
The 2D Euler equations are solved on a Sun work- 
station in double precision using a conventional finite- 
volume upwind flux-vector-splilting scheme. A C-mesh 
computational grid is used with dimensions 129x33 
grid points. High-quality lift-corrected boundary condi- 
tions are used at the far-field boundary, which is placed 
approximately five chord-lengths from the surface of 
the airfoil. 

In the present example, derivatives of three aerody- 
namic output functions are considered: C L , C D , and Cm 
( i.e., coefficients of lift, wave-drag, and pitching- 
moment, respectively). The computed steady-state val- 
ues of these aerodynamic force coefficients are given in 
Table 2. 


Table 2. Aerodynamic Force Coefficients 


C L 

+0.2830659 E+00 


+0.2070493 E-01 

Cm 

-0.2876639 E-01 


In addition, derivatives with respect to three aerody- 
namic input variables are considered. They are g (a 
geometric-shape variable), a, and M w . The geometric- 
shape variable g is simply a single arbitrarily-selected 
“y” coordinate of the computational grid on the surface 
of the airfoil. 

Calculation of SO SDs by SO Method 3 requires 
that all FO SDs are calculated first using both the for- 
ward-mode (Eq. (3) and (4)) and the reverse-mode (Eq. 
(6) and (7)) approaches. The calculated FO SDs from a 
hand-differentiated incremental-iterative (HDII) im- 
plementation of these two approaches are presented in 
Table 3, where the results are seen to agree, as ex- 
pected. 
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Table 3, First-Order Sensitivity Derivatives 



hj 

Forward-Mode 

Reverse-Mode 

dC L 

d bj 

g 

+0. 1 405406E+00 

+0. 1 405406E+00 

a 

-0. 1087323E-0I 

-0. 1 087323E-0 1 


-0.4672729E-0I 

-0.4672729E-01 

dC D 

db, 

g 

+0. 1761 807E+02 

+0. 1761807E+02 

a 

+0.1 I58625E+0I 

+0.1 158625E+OI 


-0.2382580E+0I 

-0.2382580E+0I 

dC u 

d bj 

\ 8 

+0.3171492E+0I 

+0.3171492E+01 

a 

+0.5955598E+0I 

+0.5955598E+0I 


-0. 1 603002E+0 1 

-0.1603002E+0! 


Table 4. Second-Order Sensitivity Derivatives 



bj/b k 

g 

a 


d 2 C L 

d bj db k 

g 

+0.248807E+03 

+0.250402E+03 

+0.205825E+03 

a 

+0.250402E+03 

+0. 1 84277E+05 

+0.1 6085 8E+05 

M^ 

+0.205825E+03 

+0. 1 60858E+05 

+0.1 33087E+05 


g 

+0.7 1 777E+02 

+0. 1 34379E+02 

+0.101304E+02 

d C D 
dbjdb k 

a 

+0. 134379E+02 

+0.9593 10E+03 

+0.80402 1E+03 

IVL 

+0. 1 0 1 304E+02 

+0.80402 1E+03 

+0.662088E+03 

. 1 _ 

g 

-0.663590E+02 

-0.60244 1E+02 

-0.490399E+02 

d 'C M 

dbjdb k 

a 

0 60244 15+02 

-0.449 1 98E+04 

-0.386943+04 

M. 

-0.490399E+02 

-0.386943E+04 

-0.3205 12E+04 


Table 5. Relative CPU Timings - Complete SO Method 3 


Computational Procedure 

% of Total 

Nonlinear Flow, Eq. (1) and (2) 

5.4 

Forward-Mode FO SDs, Eq. (3) and (4) 

25.5 

Reverse-Mode FO SDs, Eq. (7 ) and (9) 

69.0 

SO SDs, Eq. (22) 

0.1 

Total 

100.0 
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The FO SDs presented in Table 3 have been thoroughly 
verified for accuracy through a meticulous implementa- 
tion of the method of central finite-differences, where 
agreement to six significant digits or greater is noted in 
all comparisons. 

The SO Method 3 is implemented by application (in 
the forward-mode) of ADIFOR to appropriate pieces of 
the FORTRAN code used earlier for hand-differentiated 
forward-mode calculation of the FO SDs. The calcu- 
lated SO SDs from this implementation of SO 
Method 3 are presented in Table 4. The SO SDs of 
Table 4 have been thoroughly verified for accuracy 
through a meticulous application of central finite- 
differences applied to FO SDs obtained by the hand- 
differentiated methods previously described. Agree- 
ment to five significant digits or better is noted in this 
verification study for all SO SDs reported in 
Table 4. This verification study was not conducted us- 
ing finite-differences applied to the original nonlinear 
flow code; that approach has been documented to be 
vulnerable to severe numerical inaccuracy when SO 
SDs arc calculated 3 . The symmetry of the calculated 
SO SD shown in Table 4 is expected and results from 
the computations performed; i.e., no derivative symme- 
try was explicitly imposed on the problem. 

For the present airfoil example problem, Table 5 il- 
lustrates (in terms of percentages of the total) the 
breakdown of relative CPU timings for the important 
steps of SO Method 3 for the SO SDs. (Not included in 
Table 5 is the CPU time for the grid generation and the 
grid-sensitivity derivatives.) Table 5 illustrates clearly 
the computational efficiency of the SO Method 3 for 
SO SDs.. Recall that results of the present example are 
for three aerodynamic output functions and three input 
(design) variables, where the computational work of the 
forward-mode and reverse-mode procedures for FO 
SDs should be approximately equal (in theory, for 
hand-differentiated code, as used here). In this example, 
however. Table 5 reveals that the reverse-mode was 
much more costly than the forward-mode; apparently 
the three linear systems for the reverse-mode are stiffer 
than the three for the forward-mode. As expected, Ta- 
ble 5 shows that using an ADTFOR-assisted second 
differentiation, SO SD can be obtained extremely fast, 
if one already has both the forward- mode and reverse- 
mode FO SD. 

4.0 Conclusions 

An efficient incremental-iterative approach for dif- 
ferentiating advanced CFD flow codes has been suc- 
cessfully demonstrated on a 2D inviscid model prob- 
lem. The method employs the reverse-mode capability 
of the automatic-differentiation software tool ADIFOR 
3.0, and has been shown to yield accurate first-order 


aerodynamic sensitivity derivatives. A substantial re- 
duction in CPU lime and computer memory has been 
demonstrated by comparison with results from a 
straight-forward, black-box reverse-mode application of 
ADIFOR 3.0 to the same flow code. 

A computationally efficient ADIFOR-assisted pro- 
cedure for accurate second-order aerodynamic sensitiv- 
ity derivatives has been successfully verified on an in- 
viscid transonic lifting airfoil example problem. Accu- 
rate second derivatives (i.e., the complete Hessian ma- 
trices) of lift, wave-drag, and pitching-moment coeffi- 
cients with respect to geometric-shape, angle-of-attack, 
and freestream Mach number have been calculated. 
Second-order derivatives arc now computationally fea- 
sible, at least in 2D. 

This second-order method requires that first-order 
derivatives be calculated using both the forward (direct) 
and reverse (adjoint) procedures; then, second-order 
derivatives can be obtained in a non-iterative calcula- 
tion that is computationally very efficient. An ADIFOR 
differentiation is used to generate a number of required 
second-order terms (see Appendix) in this non-iterative 
calculation. If one already has either forward (NDV 
solutions) or reverse (NOF solutions) FO SDs, then 
upon obtaining the other FO SDs (NOF or NDV addi- 
tional solutions, respectively), one calculates all of the 
SO SDs (NOF x NDV 2 derivatives) very efficiently. 

5.0 Future Direction 

The improvement in computational efficiency 
achieved to date is substantial when the reverse-mode 
application of ADIFOR 3.0 in incremental-iterative 
form is compared with the black-box approach. Fur- 
thermore, work is presently in-progress where the tim- 
ing result for the ADII-AV scheme is projected to im- 
prove by an additional 30 percent over that reported 
herein. Thus the relative timing given in Table I for 
ADII-AV / HDII-AV is projected to drop from 4.7 to 
about 3.3, or better (close to the expected AD limit). 
This educated projection comes from a proposed strat- 
egy where the forward-pass execution of the ADIFOR - 
enhanccd, reverse-mode code will be performed only 
once (instead of during each iteration) in order to create 
the required ADIFOR log files. Thereafter, by repeat- 
edly reusing these fixed log files, only reverse-passes 
will be repeatedly executed during the iterative solution 
process for all aerodynamic output functions of interest. 

Extracting accurate sensitivity derivatives from ad- 
vanced flow codes is a challenging, computationally 
intensive task, particularly in 3D, even for the first- 
order derivatives. With the use of SO Method 3 and 
ADIFOR, however, accurate second-order derivatives 
are now computationally feasible, at least in 2D; per- 
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haps even 3D is within reach through future “in- 
parallel” implementation. 
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8,0 Appendix 

In this appendix, the terms DF' / Db k , 

DR' tj / Db k , and DG im / Db k are expanded using the 
index notation established previously. The expansion of 

DF'/Db, is 


2EL = ^Ln' + <*' +iS 

Db t 30. " * ax. ** 3b, 


3 Fj o' Q' i d Fj x' O' t ^ — O' 

— \l m y n U ^ A pj^nk T -m 


30.30, 

3 2 F 


3X„30, 

3 ! F. 


3*. 30, 


3Q.3X, 

■<LK +^rK, x ’* + 

1L- 

db k db l 


dX q dX /t 


i-<L+J^-xL + 


db t dX 


PJ 


J?IL X ' 

qk 

(24) 


In Eq. (24), the indices /,/, and A', are free, and repeated 
indices «, m, p, and g are summed. The terms of 

DR;,/Db k are obtained from Eq. (24) by replacing 
everywhere F!j with Ry and F t with Rj (and thus / 

replaces / as a free index in the resulting expressions). 


Finally, 

the 

expansion 

of 

the 1 

DG jn / Db k is 




DG„ 

9G. 

, c)G, 


3 G 

101 

7777 

a, + — ^ 

x' k 

j //77 

Db k " 

30„ 

ax (/ 

qk 

36* 


= 4, 


+ A. 




d 2 R, 


30„30,„ 


dXdQ, 


d 2 F 

1 x' k +-—±-x:,. 


(25) 


teM. 


+ A:/ 




3M0„, 36*30,,, 


In Eq. (25), the indices ?\ /«, and A are free, and repeated 
indices /, «, and </ are summed. 
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